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ABSIRACI 


In testing tne hypothesis that there is no  mncnotone 


trend in a gamma renewal process, the use of the statistic 


Yen ИУ 
J 2J 1 
where 
J 
y = 2X 
12 1=1 1 
and 
J 
Y = > 5, 
2Ј 1=1 1 
is investigated. The mean and variance of e is  devcloped 


BEES function of J and it is shown that Y is asymptotically 
J 


шаі ас J --> о for the gamma renewal process. À 
high-speed, theoretically exact gamma pseudo-random variate 
generator is developed, tested and compared kith other known 


techniques. The generator is then used to oktain the 
ШЕ ОТТОП ОЕ У through digital computer simulation for 
J 


Stall and modcrate values of J. 
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I. INTRODUCTION 


— — —. d 


А renewal process is a model for a series of events 
which occur at random times; the times between each two 
successive events are independent and identically 


distributed positive random variables. The time between 


Pwent 1-1 and event i is denoted by X , i-1,2,..., while 


1 
Brandon varıiableS represents the time to the i-th 

3i 
event. Wen тегов пе tine of initiation of the 
КӨЕ Г of events then 

i | 

ONE A (u ЛА»). 
1 n-1 n 


An equivalent representation of a series of events is 
the counting function  N(t) which represents the number of 
events occurring in the interval (0,t]. The equivalence of 


the two representations is expressed in the fundamental 


identity 
ІШКЕН санат аппетит S SEL Qu...) 
| n 
Шишо - P(N(t) < n) - РІ5 > tj), where the notation Pf } 


n 

represents the probability of the event uithin the brackets. 
Two important characteristics of any series of events 

are the so-called renewal function M(t), which is the 

EpEsected number of events up to time t, and its derivative 


the rate function: 


M(t) = EL N(t) J, 
LIES dq S. 
À (t) cr (t) 





A stationary renewal process is one in which N(t) has 
Stationary increments (see Cox and Lewis [6), Chapter 4, for 


details). Іп this case, the renewal function is 
M(t)= t / u, 


Во) - 1/1, where ú is the mean time Letween events. 

Ко EO o yi O rj cerijos rot events there аге папу 
alternatives to the renewal process.  Successive inter-event 
times may not be independent or may depend on the entire 
Euctory of the process, or the distributions of the times 
between events or the increments of N(t) may change in sone 
manner as the process continues. When these distrakutional 
changes involve the expected values of the inter-event times 
or the value of A (t), then a trend is said tco exist in the 
process. The trend may be on the serial number of the 
event; for example, a radio set may tend to ‘fail more 
freguently as a function of the total number of failures 
which have occurred. Alternatively, the trend пау be 
essentially independent of the number of events and reflect 
basically the elapsed time t; such a trend cn time might be 
okserved in the arrival of jobs at a computer facility over 
a 24-hour day where the job arrivai rate depends mainly on 
the time of day and not on the number of jobs submitted. A 
trend on time is usually expressed through a non-linear time 
dependence in M(t). 

ШЕспіс may be increasing, decreasing, cyclic or some 
combination of the three. The varicus statistical 
procedures designed to test whether an arbitrary series of 
EL has a trend of any type are quite specific as to the 
type of process and the trend model assumed. For a more 
complete discussion of tests for trend in series cf events 
(Сок апа Lewis [6], Cox [7] and Lewis [ 18}. 

The simplest form of renewal process is the Foisson 
process in which the inter-event times have the negative 


Sepcnential distribution. A more general process is the 





gamma renewal process in which the Х #5 are independent 
i 


gamma random variables with density function 

и Хх 

E(x) eee À x e ШЕ» Ор, 
ГКУ 
In this garameterization, À js a scale parameter while k 
characterizes the distribution and is called the shape 
parameter. The function T(k) is the complete gamma 
EuNEction cr integral | 
со kt =X 


ГК) = i X Ow dx... 


In testing for trend on serial number in the gamma 


renewal process the use of the statistic 


У = У Y 
J 2J / ТУ 
where 
5 
Y = X 
TJ 2 {| eet 
and 
J 
S 
244) I, 


E unvestigated in this paper. A test for trend using Y is 
J 


S I c d, and the distribution of Y is investigated 
| J 


through simulation for small values of J and through normal 


thecry results as J --»5 o. 





I BACKGROUND 
Cox [5] proposed the use of the statistic 


J 
B = 2 t 
1-1 1 


аз а test for trend in a Poisson process, where ОЕМ 


are the normalized times to events in a fixed observation 
interval of length T and J is the number of events occurring 
in (0,7], that is, the observed value of N(T). Bartholomew 
me; Gealt with the statistic 


and developed results on the asymptotic relative efficiency 
of tests for trend based on this statistic against several 
different trend models. Lewis [18] showed that the test 
based on m is not a consistent test for stationary Poisson 
processes against stationary gamma renewal alternatives, 
although it had been used for that purpose. 

The test based on the statistic B is a conditional test, 
the conditioning being on the observed value J of the random 
variable N(T). This is because the test is designed to 
detect trend regardless of the rate of occurrence of events. 
The test arises from the non-homogeneous Poisson process 
with rate function 


PE 
A(t) = expfa + bt} = e . 





=ч 


In this case іс 15 desired to test whether b=0, that is 
MU Ter there is a trend. It is possible to obtain the 
ООС Ест V L( A b) of the observations for this 
model and from Neyman-Pearson theory conclude that the test 


based on B is the uniformly most powerful conditional test 
Н : р=го,х 20 
0 
Le b#0, A 20 


Dor every value of T. For details, see Cox and Lewis [6], 
@raplers 2 and 3. 

For the foregoing model, the parameter À is a nuisance 
parameter since it has no effect on the trend; it merely 
controls the base rate of the process. It is eliminated in 
ED test by conditioning on J = N(T), this being a 
SieEricient statistic for A for any fixed value of b. 

There is an analogous situation when testing fcr trend 
on serial number and the process is observed for а fixed 


number of intervals J. Here, the trend model is one of 


exponentially distributed intervals with 


bi 
E X ] = exp{a-t bi} = е A TTT e 
1 
EI "sufficient statistic for } is nowS = > SO that the 
| J J 
Што now conditioned on the observed value of Y and the 
1d 
Extostic analogous to B is 
5 5 
Y == 2S = Г Xe- 
23 1=1 1 Te] ) 1 
V imptotically, thore 15 no difference between 


conditioning on a fixed time interval and conditioning ona 
fixed number of events but the latter turns out to be more 


convenient for the gamma renewal process. In this case 


10 





performing a test for trend analogous to the Poisson process 


test based on В requires the determination of the 


conditional distribution of M given SE 


Be THE GAMMA RENEWAL PROCESS 


In the gamma renewal process, it is both analytically 
МИН оошрешке она ту simpler to test for trend on serial 
nunber over a fixed number of observations than to test for 
trend in time over a fixed time period. The computational 
advantage of the former test arises from the ease of 
computer simulation of a fixed number of events as compared 
to the difficulties inherent in continuous time simulation 
control. Analytically, the gamma renewal process has 
certain theoretical properties which simplify consideration 
of the serial trend test. 

For the gamma renewal process with shape parameter Xx and 
the assumed trend model (Cox |71) 


bi 
ера ое (ns 25 е) 
SL 


ENS possikle to set up the likelihood function L(k,b,2AX). 
Explicit results, however, cannot be obtained for estimating 
b or testing whether b=0; it is therefore usual to use the 
Poisson test based on В іп this case. The test will 
probably be relatively powerful for k near one but not when 
К 15 small as it uses very little information about the 
individual en me C eetet iS based only on the 
centroid of the times to events and not on times between 
events. 

Elimination of the nuisance parameter X cannot be 


accomplished for the gamma renewal process by conditioning 


11 





ШІ S dq coubdnthe Poisson Case, since S is not a 
J J 


Eurticuentostatrstre for" X for all values of b. However, 


since X is a pure scale factor the ratio statistic 


J 23 19 


ЕЕ СЕ ) тог апу 5 апа К апа 15 thus a reasonable 
alternative statistic on which to base a test which is free 
me the scale factor. 

The distributions OÉ the two test statistics 
(conditicnal and ratio) are equivalent in the gamma case 
under the null hypothesis (b=0) because of a result which 
characterizes the gamma distribution.  Laha [16] and Lukacs 
[21] have established that, for a Sequence of independent 
and identically distributed random variables o SX. wd ТЕ 

n 


Mie statistics 


J 
Weed) 2% 


Ima 
— 


1 


апа 


1=1 1 


are inderendent if and only if the X #5 have the gamma 


i 
ВЕЌЕ Сга бцъјоп. Thus; iS9nuependent of Y for the gamma 
J 13 
memewWwal process. Using this result, 
БЕСІ = a учу - > 
J 23 13 i 
= Pie /Y Ss; y = у) 
22 12 19 
= р{ү СОУ = ү}, 
29 15 





КОО ОО ОК ОО т Опа and ratio Statistics have the same 
distribution. Furthermore, the gamma renewal process is the 
оу ОП он 16:5 15 the case. 

ПОШ ШО ШШ LEO: the above characterization and the 
roba bulity relationship is that, if the joint distribution 


or Y and > is asymptotically bivariate normal (a result 
13 J 


which is proved in the fcllowing Section), then д must also 


be asymptotically normal. This follows since the 


Eoudstional distribution of X given x will be 
24 24 
asymptotically normal under these conditions. Thus, 


AMO gh the asymptotic distribution of a ratio statistic 
such as = can be very difficult to obtain in general, an 
expression can be found for the gamma case. 

Trend tests using = for other types of renewal process 


could and should be considered but the gamma process was 
chosen deliberately to take advantage of the above 
Simplifying distributional results. The gamma process would 
probably not be the simplest case to consider under the 
trend alternative, however. 


Ji the succeeding sections, results on the asymptotic 
distribution of Y and the mean and variance of Y for small 
| Ј Ј 


values of J are developed. In addition, a description of a 


computer simulation experiment vhich was performed to obtain 


NS ER bution Of У for small values of J is described 
J 


and its results presented. The power of the Y test for the 
Jj 


gamma process has not been investigated nor has its power 


lative tO nonparametric tests or tests using data 


13 





— 


transformations been examined. 


NASA E IO LEAD ES TR IBUTIONS 


The statistics Y andi i for a reneual process are 
1J 24 
defined by 
у 
Y = X 
ii i=1 
Y = 5 = и] X 
AS = i' 


where oir ... X are the process inter-event times and 
J 


2 
Emnosassumed to have finite mean u and finite variance o . 


Obviously, a and > are dependent random variables 
J J 


With respective means and variances 


J 
= > Var[X | 
1-1 1 


J Var[ X і 


14 





|| 
С. 
e 
ғ 


IHM 


= Е. 


СЕЕ mj 
1-1 J: 


ML) EA 
T 


|| 
име 


1. 


име 


i | 


J 
= ENYA (Jt1 = 
а 1 


= 


= J(J+1) u. 
444210. 


2 
S = VariY 
2 | 25) 


ИРИС 


= Var. A у У 
1-1 m 
J m N 
ОЕТ) 
1= 1 3. 


= Var[X] 


к 
IMAS 
E 
к 


IC, 
{+ 
(55 
F 
lo 

| 

| 


Е Ze 


ого the correlation coefficient for the statistics Ee 


and - it is necessary to determine the first joint moment: 
J | 


m TU E s 
10% Í 1J 237 


J J 
ee cH eor CH oU M QNI 
і-і 1 1-1 1 


J 
J+ 1-k) X X 
T i K i 





J 2 
m = Ep > (IFA X. 
12 1-1 1 
E 5 (2J+2-i-k) X X J 
+ +2-1- 
i=1 k= Les i k 
J 2 
СШ C i-i) кру 
i= al 
75 $ (2J«2-i-k) PF[X X ] 
t т £ 
1=1 Е т К 
21053 ep 2. Š (isk) 
= ЕХ L * : J. T+ - 
[X ] із1 1=1 к= +1 
After some simplification this becomes 
P 2 2 
m = J (Jt FIX l+ J(J - U. 
43 449542 [e Ji 1 


meen the covariance 15 


u = Е X => NR Y = til 
11 DoE 1) (7 2) 1 


= m n 
12 1472 


2 
= 54454). С 


EN the correlation coefficient is 


= Uu ⁄ ss 
i 1 1 СЕА 


ко 


ШОШО а asymptotic results for Y as J -—> =< 
J 


necessary to shov that the joint distribution of Y 


mem asymptotically а bivariate normal distribution. 


16 


2 2 
J (Jt 1 иле. (ОТ (201 
444241. G Жақ JAEN LIE c 


z iC IS 


and Y 


24 
то do 





ШЕР 11 715 песессвату to apply the Central Limit Theorem to 


the two-dimensional vector random variable 


ы. а” 
5 VUE 


Wald and Wolfowitz [29] have shown that e is asymptotically 


normal if and only if 


is also asymptotically normal for every choice of > and >: 
To apply the Central Limit Theorem to L , define 
d 
им Е 
Jd J (Dy 


4 J 
A A A AE) Xx = EL ) 
= 1=1 1 J 


1 і=1 i 2 
J | 

= 201, +01. (9+1-1) 1 (х. 7 u) 
J 

= аи ш 


оо). 
1 1 2 
III Mate E, e 
a ЕСІ 
2 Var 
S = Var[ L 
1 2) 
= Var[ L'] 
4 


J 
I I| = "је а 
i е 


17 





Then L will be asymptotically normal if and only if the 
J 
Lindeberg condition is satisfied for every e > 0 (Cramer 


[8]: 


+ 11 1 5 Жы СЕ. y) 0 
= LM V V = • 
) 2--> о S 1=1 Bu i 


Мет Making the substitution v = a Z the integral becomes 
di 


2 2 
v dF (v) ZEE Са а је 
| а 


Ivi^es 
L 


! | 


This can be simplified by first noting that 


а . 
Пала ев” 1 
1. T 


а 72) раф тя SN 
- 1 LSR ү 


Pur Nm 97, 
1 


“p IX < ztu) 
L 


Г (2+0), 


where F(x) is the common distribution function of the X 's. 
1. 


Secondly, if 
la z[ > e s 
B 1 
then 
poe s / ја | 
1 T 
рез / Ма + и 


1 


18 





> x 44! + + J+1-1 |) 
E OA CI |! ! zi ) 


Ec o _ шу; 
1 1 2 


We can thus assume, without loss of generality, that Er 


and dida Then, 


2 2 | 2 
[ С ез: 5 2 f * nei i ШЕН t. 
EV 1 l zi>es + 
ее. | 2 | a > 
J 
> v)} 
E RA CM n 
J 2 2 
ор. ва | | ZIGE (ZEN) 
55 1211 М И +i J) 
1 1 T? 
2 
02-2) SLM ОЕ ау. 
G Z ` С. 


Since by hypothesis the X 's have a finite second moment and 


це 


since 

lim es Zi UIS 

Jo-> œ 1 ( Ine | 
> (е/у ү 11ї (5 /4) 
| 2 Ј--> с 1 

Ж” 21- 
2 ес lim 1 > а 
( А 2! 2-->о 9 lie, d 
— ( 1 IET š 1 1 +1] 1-3 Wie 
> e im + Ј+1-1 
c J 5 ГЕС C£. m [ I Е ju 

> co 


the integral (2-2) approaches zero and so the  Lindeberg 


condition (2-1) is satisfied. 


T9 





It has thus been shown that Y and ү аге 
1Ј 2Ј 


asymptotically jointly normal. From the known properties of 


the two-dimensional normal distribution (Gnedenko [10]), the 


conditional distribution of Es given E nust also be 
J C 


asymptotically normal with mean and variance 


| 
= 
+ 
о 
Y 
~ 
tn 
—. 
ма 
| 
= 


FLY [Y = 
EY) 2 2 1 1 


D 2 
Vari y Y = = 5 jo 
[ I Ж Ya PE o ) 


2 
4 J tse ty Case ty о 


= Jd(d -1) c 
t 2 


ШОШ the results of the previous Section, 22 15 а15о 


asymptotically normal with mean and variance 


EY = El Y Y Ү = 
[ 21 24 2 14 | 14 Y 1 
= 
2 
Varl r =. Уаг[ У ү y = 
[ 2] [ 24 2 ПА | 2Ј Y 1 


J d 1 i 
== -— С с 
did 54) ( / У.) 


The foregoing result then suggests the following 


procedure to test for trend in a gamma renewal process. The 


20 





observed values of ^ and © can be used to accept or 
J Ј 


Be ect the mull hypothesis of a trend-less process at any 
desired significance level a using the conditional normal 
distribution given above. Rejection of the no-trend 


Iypothesis will take place if the observed value of Y is 
J 


Er ter than the 1-« /2 quantile of the distribution of 2 


(eer increasing trend) or if i is less than the «/2 
quantile (decreasing trend). 

For the distribution of the sum of a series of 
independent exponential random variables CO become 
approximately normal, however, requires over 20 terms; the 
gamma renewal process with shape parameter less than one 
requires several terms just to produce a single exponential 
variable. Thus, on the order of 100 events could be 
required to use the asymptotic test. For a small sample 
(roughly, one with fewer than 207k events) it may not be 
possible to use a normal approximation to the distribution 
of the test statistic. 

If sample size is so small that the asymptotic result 


cannot be used, it is then necessary to develop the 
ПИ оосјог о: У in order to use this test for trend. If 
J 


Poe distribution cannot be determined analytically, computer 
Simulation must be resorted to; it is still possible, 


however, to use the probability relationship between the 


w tribution of Y and Y given Y to obtain the moments 
J J 13 


"LY. 
J 


21 





ПИ p Ir P АМО VARIANCE OF Y FOR SMALL SAMPLES 
Ј 


Fron the foregoing, it can be concluded that in a gamma 


renewal process 


II 
[т] 
сан 
"d 
NS 
к< 
Li 


d 


uU / o 


 іпсе Үү /Y is independent of Y . Thus, 
2da lJ 19 


PI Y = m 
[ m / 
zu 
2 
лаг] у, | 
ЕГУ 2. = gr v^ / y 7 
J EN Їз 


2 2 
uns E 

2 2 2 2 
C има Sk + E 


II 


2 2 
dtu} 2 (2J+ 1) G + 3J (J+1)u . 
ap e 


Introducing the mean and variance of the 


r able, 


с 
II 


ких 


2 2 


2 2 
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gamma random 





2 
G = К ИЛ 


2 2 
Е[Ү ] = Jet 202541) ких ух 32 941) (К E: 
poU А е; Ef 
И ОК. 
УШ 202941 атата) 


u rC inal result for the variance is then 


2 2 
Var[Y ] E[ Y] а се 


Je ee 
HEUS ЕТ 


Note that for the Poisson process (for which k=1.0), 
Var[Y ] - 4,1 blocca wwecsresublt' for the sum of 4-1 
d . 


uniform random variables. When k#1.0 then, approximately, 


~ 
c 


о 4 C (xy, 
( 2 T ( X) 


where C (X) is the coefficient of variation (which is 1/k 


for the gamma random variable.) Thus, -the greater 


variability in the inter-event times is reflected directly 


Mere increased variability of the statistic RE rrom the 


point of view of a trend test, the variability of ic is 


increased to allow for the possible confusion of an actual 
trend with a sequence of long intervals in the stationary 
process. As mentioned previously, this indicates that more 
powerful tests for trend could be developed which use more 
information about the interval distribution than just the 
coefficient of variation. 


The succeeding sections discuss determination ОҒ the 
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Ug stribution of ze for small samples through digital 


computer simulation. Values of the predicted mean and 
variance of E were calculated for each simulation run ana 
the results are tabulated as follows (note that the mean 


does not depend on k): 


J MEAN VARIANCE 
КЕСУІ Т” 50 К = .75 К = 1.0 


10 DD 4.125 22375 125 75 I 5/45) 
80. 15.5 18.5229 8.814 4.682 3.188 2 17 
8025.5 34.708 15.428 8.010 5.409 4.083 
lU 50.5 15. 150 32.048 16.338 10.964 2222) 
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ieee ои рака GANMAZVARTATES FOR SIMULATION 
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Although there were several known methods of generating 
gamma random variates (see Section ІТІ.Е below), all were 
judged deficient for large-scale computer simulation use; 
some are statistically inadequate for large samples and all 
require too much computer time. A new method was therefore 


developed to overcome these difficulties. 


ШІ DESCRIPTION OF ALGORITHM 


Marsaglia's rectangle-wedge-tail methods for generating 
normal [ 14,23} and exponential [22] random variates were 
modified to produce unit gamma variates with shape parameter 


K less than one. The basic approach is to select with 


respective probabilities ЖЫ, "и one oft а series of 
n 


Andon variables with distribution functions тылы К 


ENDE (x). If x is the generated random number then 
n 


E E AC AA + P F (x). 
1 1 ДА | n a" 
Taking derivatives, 
eo) werner... ЕР f (x). 
) 15,6) + 2,5, 00 Е 09) 


The method is then equivalent to a geometric decomposition 

EE Uc density function of the variate to be generated. For 

 Ісіепсү, ілер 's and F 's are chosen so that most of the 
at ШІ 


time the generated variate 15 selected from a simple 
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distribution (preferably a uniform distribution). 

Gamma variates with any desired scale parameter are 
obtained by multiplying a unit gamma variate (that is, one 
with scale parameter one) by the reciprocal of the desired 
Seale parameter А. Thus, the method developed generates 
only unit gamma variates, which have density function 

k-i =x 
Ғ(ху- 1 X e š 
ПЕТ 
The unit variates can then be scaled to produce the desired 
result. 

It is difficult to apply the decomposition technique to 

the gamma distribution because the density function has a 


singularity at the origin when the shape parameter is less 
than one. А low value of x (called m is thus selected as 
the starting point for the decomposition; values of the 
gamma variate less than 22 are generated with probability P 


by means of a series approximation to the inverse of the 


incomplete gamma function (see Subsection  III.B. 1.) À 


сс ог x 's is then determined by the relations 


i 
но = E 
1 1 
(3-1) ћ - 2h. 
1 = 
X = х, ou Ее. Г). 
1 1 Í 1-1 


This procedure continues until the area under the density 
mmetlOn Curve from zero to x is greater than 0.999; the 
i 


value N is then defined to be 1-1. Dependent upon k, the 


parameter N varies between 20 and 100. Values of the 
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variate greater than x are generated by an iterative 
N+ 


psocedure (Subsection III.B.2) with probability p 1° 
+ 


The density function has thus been decomposed into N-«2 


regions: a "head" region below B. а "tail" region above 


X and a series of N vertical strips. Each strip is 
N*1 


divided into a large rectangle and a Small wedge, as shown 


in Figure 1. fhe probability for the strip as a whole can 


Е (х) | 





s 
— == аза 2. 


X. X, 
1 1%1 


| 
| 
Í 
v EMT Soc ER: 
— 
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Figure i: Decomposition of strip i. 
be found by using the incomplete gamma function, 


СУ ES y 


g (K; x) y e dy, 


= 1 | 
ШЕТ 0 
ШЕ Ол їз readily evaluated by an infinite series [1] 


(К: х) 1 5 Wc 
JUN A] — usse 2 = 
PAK) n=0 iT у- 

suitably truncated for computational purposes. 


ШЕБЕР” 15 спе probability for the strip as a whole, p. 
1 1 


the probability of ihe rectangle aná q the probability of 
i 
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ES) 


the vedge then the following relations exist: 


ах уе guo er 


1 1+1 
pee EX ), 
ДЕ 1 1%1 
а == Tp 
1 1 1 


Поште values o, q ор because of the extreme 
1. 1 


steepness of f(x) near the origin; this accounts for the 


` 


 Есасіпп values of h in (3-1). For i greater than 5, 
1 


this relationship reverses so that most of the area under 
the f(x) curve lies within the rectangles. In fact, it vas 
found empirically that 


N 5 
> p. = 0.172 


121 1 
Ber most values of к in the range 0.1 to 1.0. Thus, the 
gamma random variable can be sampled from а uniform 
distribution over 70 per cent of the time. Presumably, an 


increase in this probability would result from choosing each 


h optimally, but, the gains to be made are small. 
1 


1. Binary Search Schene 
With the decomposition complete, there remain two 
problems: efficient selection of one of some 100 sampling 
distributions and methods for sampling from each of the four 
meererent types “of distribution (head, tail, rectangle and 
wedge). 
To formalize the selection problem, suppose that 


events Eo ... ‚Е are to be selected with probabilities 
M 
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P sP r œ.œ’ ‚р ‚ respectively. In the gamma generator, М = 
М 


Ра. A cumulative probability vector P = te ay ees Р) 


is defined by 


i 
P = Ёр. ME 272. + М) • 
1. J=1 j - 

ima candcm number U uniformly distributed on (0,1) is 


generated, the selection of an event can be carried out by 


comparing U serially with P wand stopping when, for 


the first time, U < P for some n; event n is then selected. 
n 


This methcd requires on the order of  M/2 comparisons and 
would be far too slow to solve the selection problem. 
Substantial saving of comparisons results from 
applying the binary search methcd to the problem of finding 
Such that P COE T pus For this scheme, an interval of 
= | n 


uncertainty is defined by selection of the indices i and j 


EE 0 SU <P “at all times; initially, 1 is taken to 
1 2 


be zero (with E - 0) and j to be M. U is then compared 


with and the interval of uncertainty adjusted 


BD 
[ (1*3) 72] 
 Әтіппа to the result. This procedure is continued until 


BRSEdifterence i-j is equal to 1; then n= j. It is not 
difficult tc show that с! comparisons are needed to find 


n, So that this method 15 considerably faster than the 
sequential search described previously. 


II saving Ot Comparisons is still possible, 


meer eS L nce the exact distribution of the p 's is known, 
i 
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the search procedure can be 'skewed" so that more probable 
events can be found in fewer comparisons, thus minimizing 


the average number of comparisons. This can be done Бу 


прагіпа U with the value of E which is closest to tne 


average of P and P and proceeding as in binary search. 
1 2 

Finding the value of К for given i and j is no longer as 

Per however; in the actual method used, two arrays of 


inđices called Last and Next are precalculated. If U is 


less than P the following comparison is made with P 5 
k Last (k) 
if U is greater than P then P is used for the 
k Next (k) 
succeeding comparison. Termination of the search is 
indicated when Last(k) = 0 (choose event k) or Next(k) - 0 


(choose event k+1). 

The binary search scheme used here is closely 
related to Huffman (minimum length) codes in infcrmation 
thecry and is essentially similar to Knuth's optimum binary 
БЕСІСП tree algorithm [{ 15]. а Бе shown [15] that 
between R(E) and H(P)*2 comparisons are required for this 


method, vhere H(P) is Shannon's entropy function 


H 
H(P) - .2 P. log (1/p.) 
1= 2 2 1 
mers example, for k=0.1 it was found that M = 2N + 2 = 164 


and that less than 6.4 comparisons were needed on the 


average to select an event. 
2. Rectangle Method 
If the binary search scheme selects rectangle i and 


if Vis a second uniform (0,1) random variable independent 
OO, then 
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is a random number selected from the i-th rectangle. 


3. Rejection Method for Wedges 
О она а [23] discuss in some 
detail a method for sampling a random variable whose density 
munction is almost linear. This method is used without 
modification to sample from the wedges. 
mes talves ok a апа П were found for each wedge 
Д 1 
ae the gamma density function f(x) lies within the 


Ap indjeated in Figure 2, that is so that 


CeO Sean eee = s p) о (Xx х )/h. 
1 1 1 i 1 a 1. 1 


тон о = у ру, 
1. | 


Smee f(x) is concave upwards for all x and since the 





м 


a ШӘ ты 
i 


Figure 2. Sampling from the i-th wedge. 


ВО ћи is fastest for a, close to b , the lover line is 
o. 1 


taken to be the tangent to the density function curve 
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mPemallel to the chord from (X,  ,y.) ao AL. The 
Iti 1 л т 


tangent can be found by Newton-Raphson iteration. 

The wedge algorithm is then: 

(1) Generate tuo uniform random numbers, U and V; 
if С > У then exchange U and V. 


(S (C 2 to x + Un + 
I j^ 


БИЕНІ OS E = (а ту J/(b =y ] go to step (5). 
i J т de 42 


(4) Tf V > U + £(4)7b go back to step (1). 
1 
(5) 2 15 the desired variate; stop. 
Hhen a and b are close together then r will be 
i al i 
nearly опе and the algorithm will terminate as a result of 
Step (3) mucst of the time. Fer most of the wedges in the 
Mamma generator, values of r are greater than .65 sc that 
i 
NA) needs tc he evaluated (step (4) ) less than 1/3 of the 
time. 


Knuth | 123] gives a proof that this method works 
porerly. 


Pee APERCAINATIONS FOR THE INVERSE GAMMA CDF 


Any continuous random variable can Ee sampled if its 


Merse distribution function F can be found, for if U is 


a uniform (0,1) variate then Z = F (U) is a variate fron 
the desired distribution. Unfortunately, there is no simple 
inverse for ‘the gamma distribution function with shape 
Parameter less than one; Phillips and Beightler [26] 


attempted a rational approximation for the inverse vith 
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mettle success (see Subsection  III.E.2). In the present 
method, two different inverse approximations are used for 
the head and tail regions. Но inverse is needed elsewhere 
Since the wedge sampling technique described above uses 


f(e), the gamma density function, and not the gamma inverse. 


e Low Values of Z 


a. Method 


The series expression for the inconplete ganna 


function (gamma CDF) can be expanded as follows: 


Ba) = g (KX) 
EE F(x) = 1 5 x 
E ° (X) = = X 
EE S) Ary 
k 2 
cx cmq нв и Xe ee Je 
ТЕСЕ K Jer КЕНТ 
If Z = F(x) then 
k 2 
Kh e x 1 = Ke KE ANNA: 255 
k*1 Z(k*2) 
1/k 2 1/k 
fee (k+1) 2 | А uctus td / р 


ES E Ку X 
K+T РЕТИ 
When AU << 1 the first-order approximation to 


enhe inverse, Т сап be used. “This corresponds to the case 


when the variable 15 to be sanpled from that part cf the 


head regicn close to zero. 


1/k 
['(k+1)z] 2 \ 


> 
| 


сет! 
P (2) 


E 





When the E X term cannot be ignored, a better 
+ 


approximation is obtained by solving 


1 k O a 
Xx = X EN NE x 
1 i k+T | 


Ы 2 
м ко olx 
: K К+Т E) 


= 2 | 3 
- x — T x + O(x ) 
Кет 


for the second-order approximation, B The result is 


uc OE E E ] 


-1 
Frl. 


À more convenient form of 2 for computation is 


uns AL A 1) ]. 


DEE aC y 


The value of E (probability of the head 


meogron) is determined empirically for the gamma generator by 


the formula 


во that the accuracy of the approximation to the inverse 


incomplete gamma function vas investigated for values of 7 
in the range zero to р for various values of k. It was 
0 


found that the approximation improved as 2 --> 0 and as 
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k --> 0: in no case did the relative error exceed 10 in 


double precision or 10 іп single precision. 


28 ШЕ ог 2 Near One 


a. Method 


Newton-Raphson iteration is used for 
approximating the inverse gamma distribution function for 


values of Z near one; that is the equation 
век = 2 = 0 


is solved for x by using Neuton's formula 


x ЗА) a 
п +1 П П n 


Six-digit convergence can be acheived after four or five 
Meerations Of this method starting with E = 1.0. 


Application of the method requires evaluation 
БЕ (х) апа h'(x). The latter is just the gamma density 


СЕ оп 


hec mix = 1 р = 
г (К) 
while the former can be found by summing the series (3-2). 
Since x is likely to be very large when Z is near one, 
however, as many as 30 terms of the series must be summed 
for convergence with a resultant loss of speed and the 
cumulaticn of serious round-off errors. For this reason, 
g(k;x) is evaluated by a continued fraction approximation 
[1] vhich has the added advantage of producing the value of 
ШООСО(К: х). Since the method is applied only to values of 2 
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"T". 


in the range (0.999,1.000) the ability to deal with 1-2 


[eras directly to a net cain in significant digits. 
БӘР Accuracy 


пита u I (T но amvestigqate the accuracy of 
this approximation in the same way as for the previous 
Eod because of the loss in significant digits in 1 - Z as 
2 --> 1. A series of values of y in the range .00001 to 
ШО was defined; for each value the approximaticn was 

> 1 Ë 

E to finda x= F (1-у) апа then z = F(x) was obtained by 
ne continued fraction approximation. The relative error 


between z and 1-y vas in most cases Zero and in no case 


greater than Heí10 . 


DESCRIPTION OF COMPUTER PROGRAM 


The gamma generator ма 5 implemented as a 
EU UN-callable subroutine with each call returning a 
single floating point (type REAL) gamma variate. One 
version was written in the IBM implementation of FORTRAN IV 
while another vas written in a combination of FORTRAN and 
EM 560/Assembler F. Copies of listings for both prcgrams 
are on pages 77-84 and 85-88. 

The basic approach in both generators is to divide 
the subroutine into two parts: a relatively lengthy (and 
Slow) initialization section which calculates tables and 
E-tants aud а short, fast section which is called 
repeatedly to handle the actual generation process. The 
initialization part of the program was given the name GMINIT 


(for GaMma  INITialization) while the generator was called 
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GAMA. Calling sequences and arguments are discussed in the 
following Subsection. 

GMINIT accepts values of k in the range 0.05 sk < 
200. Tt produces the cumulative vector PROB (analogous to 
Вт Subsection III.A.1); before being placed into PROE the 
component probabilities are ordered with the smallest first 
to minimize round-off errors.  GMINIT also finds the arrays 
NEXT and LAST for use in the binary search method as well as 
ee arrays %,H,R and B for sampling from the rectangles and 
vedges. Finally, GMINIT calculates constants for the head 
and tail approximations. Sample output from GMINIT is on 
pages 70-72. 

Two geurther subroutines Were also written: 
IGAM(K,X) which evaluates the incomplete gamma function 
g(k;x) and INVGAM (K,Z) which uses Newton-Raphson iteration 
to solve the equation g(k;x) = 1-2 for x. INVGAM was also 
programmed in both FORTRAN and Assembly language. Listings 
meee both ІСЛАМ and INVGAM are on pages 89-90 and 91-95, 
respectively. 

In addition, a uniform random number generator is 
required. Since an average of about 2.5 uniform numbers are 
required for each gamma variate, timing characteristics of 
the uniform generator are a significant factor in the 
performance of the gamma generator. The particular 
generator used is called RANDOM and was developed by Lewis, 
Goodman and Miller [19]. RANDOM is very fast; it is capable 
of generating a uniform deviate in about 15 microseccnds. 

The generating routine, GAMA, carries out binary 
search with a Single uniform random number and then samples 
from the proper sub-distribution using one or more 
additional independent uniform variates. It also performs 
scaling on the generated variate by multiplying it by a 
scale factor BETA which is an input parameter to GMINIT. 

For small values of k (less than 0.1) there is some 
Beorabılıty of a fixed point underflow error, i.e. that the 


generated  variate will be less than the smallest possible 


5) 





— 


СО 
floating point number (about 10 Cor he ВИ 360). Тһе 


FORTRAN standard fixup for this error is to return a resuit 
of zero, which might be satisfactory for some applications; 
it could create serious problems, however, if, for example, 
a log transformation of the gamma variates is required. If 
a large enough scale factor is used, modification of GAMA to 
per form prescaling before applying the low value 


approximation would reduce the probiem. 


2. Implementations 


a. FORTRAN IV Version 


In the FORTRAN version of the generator the 


instruction sequence 


CALL GMINIT(K¿BETA) 


CALL GAMA(Z,1X) 


results in the initialization of the generator for shape 
parameter К and scale parameter Aà = 1/BETA and the 
assignment of a gamma random value to Z. IX is an integer 
seed for the random number generator RANDOM and should not 
normally be modified in the main program. The same sequence 
of gamma deviates can be reproduced at any time, houever, by 
EL Unitrialiizing IX to its original value. 

IBH FORTRAN IV has several features that are 
not in American National Standard FORTRAN and which may not 
be implemented in other compilers. Two of these are used in 
the generator and its subroutines: the built-in function 
GAMMA(X) which computes the gamma function of its argument 
and the ENTRY statement which allows a subroutine to have 
 ҮҮСІріс entry points. The ENTRY option is used to cut down 
Oh overhead in differentiating which of the twe parts 


(GMINIT ог GANA) of the subroutine is being called; it is 


38 


о-ы | | | 





equivalent to having two independent subroutines sharing a 
COMMON block. Besides modification o these areas, 
conversion to another conputer would require adjustment of 


convergence test values in INVGAM, IGA4 and GNINIT. 
b. Assembler Version 


The Assembler version of the generator was 


written to save execution time in GAMA; the instruction 


sequence 
CRELISSTRRIEIK,BETA,IX) 
2 = САМА (0) 
has the same effect as the previous sequence. Linkage 


overhead in САМА has been substantially reduced, however, 
since no arguments need be passed (the zero in the 
instruction above is a dummy argument and is required by 
BESURTRAN conventions). The capability of re-initializing the 
seed for RANDOM was not retained in this version. 

START is a dummy section of GAMA which calls on 


modified version of GHINIT. Since execution time was a 


secondary consideration for the initialization routine 
(which was normally called only once), СМІМІТ was 
maintained in FORTRAN for ease of programming. The only 


Meditications to GMINIT in this case are in the calling 
sequence and in the values produced for the arrays NEXT and 


LAST, which were changed to allow for easier search. 


3. Timing and Core Requirements 

The generator was written with the sole aim of 
achieving the fastest possible execution time without regard 
for storage requirements. The execution times obtained were 
highly dependent on k; for example, more rectangles are 
needed for louer values of k so that the binary search 


nethod occupies more time in these cases. Samples of 10,000 





— 


variates were run on the Naval Postgraduate School IBM 
360/67 computer and the generation times recorded for 
several different values of k; the average times per variate 


(in microseconds) were as follows: 


k FORTRAN IV Version Assembler Johnk's 

G Compiler H Compiler Version Method 
0.10 430 327 201 824 
025 Зао 310 29 864 
0550 350 Zu 186 889 
Ш 75 339 213 178 SS 
re 325 270 179 803 
АУ егаде 366 291 187 851 

By way of comparison, an implementation of 


Marsaglia's normal generator on the Same computer requires 
about 60 microseconds per variate while a Marsaglia 
exponential generator takes 70 microseconds per variate. 
Nehnk*s method is discussed in Subsection III. E. 3; TE IMSS 
included here for purposes of comparison. 
| Although this method of evaluating run times 15 
easy to implenent and gives a realistic estimate of what can 
be expected in practice, it is somewhat misleading since 
competition Tor systen resources between jobs in a 
multiprogrammed computer causes considerable variability in 
observed run times for a given job. The observations above 
were recorded when the system vas relatively idle; on a busy 
afternoon these times may increase by as much as 15 per 
cent. | 

Bu ccu MEroZcombiete the set-up phase in GMINIT 
were also highly dependent on k; the following run time 


values were observed: 


k Set-up Time 
(Seconds) 
О fous 
0225 .348 
0.50 E212 
OTRO ALGO 
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Since a much larger number of rectangles is required for 
decomposition when K is small, it can be seen that the 
additional computation time for set-up in these cases can be 
Buster costly. ји fact, enough time is used to generate 
between 900 and 4500 variates, depending upon the value of 
k. 

Both the Assembler and FORTRAN versions of the 
generator require excessive core storage, at least in 
comparison with other methods for generating random variates 
(see, for example, Ahrens and Dieter [2)). At the Naval 
Postgraduate School installation, where the minimum job core 
allotment is 100,000 bytes (25,000 words), the size of the 
program was not а restriction, although it could be at 
another installation. Both versions of the subroutine 
require, in addition to space for their oun code and the 
motes, core storage for the FORTRAN built-in functions EXP, 
DEXDP, MOG, DLOG, GAMMA, GAMMA and SORT, as well as the 
user subroutines IGAM, INVGAM and RANDOM. Core requirements 


(in bytes) for the various subroutines are as follous: 


FORTRAN IV Version Assembler 
G Compiler H Compiler Version 
САМА 13206 12400 5664 
GMINIT 6860 
INVGAN 1070 866 684 
ТСАМ 740 574 740 
RANDCH | 320 320 320 
Functions 5056 5056 5056 
Total 20392 19216 19320 


Further improvements in both space requirements and 
run time may still be possible through the use of in-line 
ESdung instead of subroutine calls. These modifications 
would have the further result of making the generator 
independent of vendor supplied mathematical  subroutines. 
Ahrens and Dieter [2] point out some: advantages of this 
approach in their improvement On Marsaglia's normal 


MÉnerator. 
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ПК SIATISTITCAL TESTING 


Having thoroughly investigated the numerical accuracy of 
the various approximations used in the generator, it was 
then necessary to test the generated  variates for 
"Statistical accuracy". Formally, the null hypothesis 


H : The generated deviates gu me were 


sampled from a unit gamma distribution 


with shape parameter k. 


was to be tested against the alternative 


H : The generated deviates were not sampled 
from the unit ganma distribution, 
ОЕ the wide variety of standard qoodness-of-fit tests, two 
were chosen for testing s the chi-square test and the 


Anderson-Darling test. Inemadag)stbion ues reproducttvre 


property of the gamma distribution (that is, that the sum of 
Ено gamma random variables with shape parameters se апа К 

and equal scale parameters is also a gamma random variable 
EG Shape parameter А was exploited in a test on the 


sums of several variates. 


1.  Chizsquare Test 
A modification of the chi-square test suggested by 
Naylor [25] for uniform random number generators vas carried 
BEN Por this test, a set of n quantiles of the cbject 


EuStrabution (1.e., the gamma distribution) is first found; 
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the i-th quantile ris the solution to the equation 
1. 


п 
i 
Quantiles for the gamma distribution are unique. А sanple 
of size N variates is then generated and the number of 
Maitiates f falling into each interval г. а Ep is 


i i-1 1 


recorded. The stati- tic 


2 3 В 2 
X = n = N 
1 N 21 i n 
2 ° 
is then found for the sanple; E. has approximately a 
chi-square distribution with n-1 degrees of freedom. This 


procedure is then repeated M times and a second statistic 
2 . . - Ф “ 

X 1s calculəted in an analogous manner, using m  quantiles 

of the chi-square distribution with n-1 degrees of freedom, 


2 
The value of > is then used to accept or reject MOS 


Values chosen for the actual test vere m = n = 10, 
N = 1000 апа М = 100. Thus, each test involved generation 
ШИ 100,000 variates.  Tabulated results for various values 


of k are as foliows: 


k X ° Chi-square Minimum X : Maximum X Í 
2 Level 1 1 
"06 4.8 . 149 2.24 о 
. 10 0.6 32 1.80 222220 
515 122 ео 1.90 22205 
E50 4.4 . 117 1.66 | DTI 
290 isa SA 166 29.059 
895 gc . 544 ШІ) 21298 


43 





— 


The column headed "Chi-square Level" gives the quantile of 


the chi-square distribution with nine degrees of freedon at 
2 - 
the observed value of i . Based on these results, m could 


not be rejected. 

There are some serious drawbacks to the chi-square 
test, not the least of which is the fact that the statistic 
depends strongly on the specific quantiles taken and is not 
an invariant for a particular sample. Nevertheless, the 
test proved to be quite sensitive even to minor variations 
in the generator. For example, it readily detected the 
Statistical effect of the omission of a Single wedge (whose 
probability was .002) by the binary search procedure, 
producing composite scores greater than 50 in this case. AS 
is indicated in the sequel (Subsection III.E.2), the test 
was also used to reject the Phillips generator, which 


acheived scores greater than 800 for small values of k. 


De Anderson-Dariing Test 


a ww A A о сы A A «=> mn 


The test of goodness-of-fit proposed by Anderson 
mea Darling [3] 1S a distribution-free test in that instead 
of the observed values of the variate AA a the 

n 


values 
uU. 75% 
1 1 


are used. j 15 factethe distribution function of 


mer xX 's (that is, if the null hypothesis is true), then the 
1 


ü s are ПЕР О Size п Егоп а uniform (0,1) 
i 


IA 


E ruüubution. If the u 's are now ordered so that u 
M (1) 
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u I S UU uc Psemserson-Darliing statistic 


(2) (п) 


Е 
a (el ) in a) 


= 
T 
| 
= 
! 
1 
IMS 


n end 
ТН Ји (1 = ч ) ] 
(1) 
can be computed. The null hypothesis can then be accepted 
2 


ШЕ ОН is not too large. Levis [17] has determined the 
n 


2 
distribution of Y for various values of n. 
n 

It should be mentioned that the Anderson-Darling 
statistic gives the most weight to observations from the 
tails of the distribution, where u is very close to Zero 

(1) 

Or one. Thus, the test was ideal for the gamma generator, 
since the rectangle-wedqe technique had been successfully 
applied ina wide number of other cases and the main 
uncertainty for the generator lay in the statistical 
properties of the head and tail approximations. 

Accordingly, Anderson-Darling statistics were found 
for 100 samples of 100 variates each. The results of Lewis 


ШІ? | меге then ЕО опа ten  Guantiles of the 
NS 
distribution of e SO that a chi-square statistic could 


be computed for the sample of 100 Anderson-Darling 


Eatrstics. The results were as follows: 


2 2 
К Chi-square Chi-square Minimum W Maximum Y 
Statıstie Level n 
.10 4.269 ЗО, Т 3a 500 
50 1:295 slow 2175 5.071 
7 90 4.306 „110 22:21 2292 
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Шпсе again it was decided to accept ie on the basis of the 


test. 
3. "Reproductive" Property Test 


It is uell-knovn (Cramer [8] that the sum of two 


independent gamma random variables with identical scale 


parameters and respective shape parameters S and z also 


has a gamma distribution with shape parameter n = Since 


~l 


a gamma random variable with k=1.0 15 an exponential 
variable, a possible test for e is to investigate whether 


Me statistic 


м о, 
n 1= 1 1 


Mas а unit exponential distribution, where the z "s are 
i 


generated by GAMA vith k = 1/n. 

The test was carried out for samples of 100 V 's 

n 
МИ - 20,68, 10. For each of four such samples, an 
Anderson-Darling statistic was computed and the maximum 
value of the statistic recorded. The approximate 
significance level of the statistic was determined by linear 
interpolation in the table of the asymptotic distribution of 
2 


W given by Lewis [17]. The results are as follows: 
n 
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2 
n k Maximum W | сЕ nate 
100 Significance Level 


2 .500 1.577 2841 
4.250 2.045 .913 
6 .167 1.031 .658 
8 .125 1.469 .816 
10 .100 1.525 ..829 


As a further test of whether the statistics V have 
n 


wn exponential distribution, the statistic 


ve = Minimum (V V ава V ) 
n nod n;2' I n; 1000 


was investigated; this was done specifically to test the 


statistical effect of the low value approximation to F (Z). 


It is well-known that the minimum of m unit exponential 

random variables is an exponential random variable with mean 

1/m. Thus, samples of 100 observations of V' were generated 
n 

BOL n=2,4,8 aná an Anderson-Darling statistic computed using 

the exponential distribution function with mean .001. The 


results are as follows: 


n k W | M eate 
100 Significance Level 

2 0.500 . 903 EBD 

4 0.250 . 431 . 182 

Ce Un VEAS Y „936 


Às can be seen, Ыс is readily accepted based on these tests. 
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Bee COMPARISON WITH OTHER METHODS 


1. Naylor's Method 


The method of Naylor [25] for producing gamma 
variates with non-integral shape parameter was not seriously 
considered for the simulation since it is not defined for 
k < 1.0. The method involves sampling a mixture of two 
gamma random variables with integral shape parameters tk5 
and tkJ*1 (the notation tk! denotes the truncated integer 
part of К). Berman has shown [28] that the method is 
Satisfactory for small simulations wben k25 but his results 
demonstrate that the method is not statisticaily acceptable 
for sample sizes on the order of 100,000 for any value of k. 

The timing results for the Naylor method obtained 


by Berman [28] for an IBA 360/65 computer suggest that a 


better vay to generate non-integral gamma random numbers is 
to add an integral gamma variate with shape parameter tk. 
and a variate from GAMA with shape parameter k - “Ки. For 


EM this modification adds less than 20 per cent to the 


execution time but renders the result statistically exact. 
2. Phillips's Nethod 


Phillips and Beightler [26] proposed the use of a 
rational approximation to the inverse gamma distribution 
function for generating gamma variates; three different 
approximations were developed паса - 7205 220 < kis 5.0 
НЕНІ к > 5.0. The technique was extensively tested by 
EULUps using the Koimogorov-Smirnov goodness of fit test; 
unfortunately, this test is notoriously insensitive to 
deviations in the tails of the distribution. Тһе chi-square 


goodness of fit test was therefore applied to the method for 
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К < 1.0 and it was found that the hypothesis that the 
Phillips numbers were sampled from the gamma distribution 
could be rejected at any desired confidence level 
(chi-square scores in excess of 800 on 9 degrees of freedom 
were observed). The method was not investigated for k>1.0. 
Furthernore, this nethod reguired substantially more 
execution time than the method developed here (on the order 


of 600 microseconds per variate). 
3. | Johnk's Hethod 


Johnk's method ps was the only known 
theoretically exact means of generating gamma variates with 
non-integral shape parameter. It is a rejection technique, 
meme the wedge method of of Subsection III.A.3. Por 
generating gamma variates with shape parameter k less than 
1.0 the method can be set forth as follows: 


(1) Generate two independent uniform (0,1) random 

variates U and V. 
р I (T h) 

(2) Set K = U and Y = V s 

Е ош кану O. go. back to (1). 

(4) Generate a new uniform variate Е апа set 
IS ln R. 

(5) 2 = We X / (X*Y) is a unit gamma variate. Scale 


Tf required and stop. 


Eronslatrons of  Johnk's original proofs of the correctness 
of this methcd are contained in the Phillips and  Beightler 
[26] and Fox [9] papers. ‘Essentially the method generates a 
beta variate and then uses the Lukacs-Laha results [16,21] 
to produce a gamma variate.  Johnk shows that when k is less 
than 1.0 the rejection in step (3) will occur between 1.00 
and 1.27 times per variate on the average. 

Johnk!s method is not satisfactory for large-scale 


simulation se because it requires excessive time; see 
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Subsection III.C.3 for the observed times for a FORTRAN 
implementation of the method. The time can be improved by 
recognizing that if U is a uniform (0,1) variate then -in U 
is a unit exponential variate. The time to generate U and 
evaluate its natural logarithm is about 130 microseconds оп 
the IBM 3607/67 computer while an assembler implementation of 
Marsaglia's exponential algorithm can deliver an exponential 
meerate in about 70 microseconds. The improved algorithm is 


then 


(1) Generate two independent exponential variates, R 
апа S. 
=R/k SEIN) 
(2) Set X = e and Y = e : 
ПИО Ин OC NY ^T do back to (1). 
(4) Generate a new exponential variate T. 
Coy Deliver Z =T*% ° X /(XtY). Stop. 


Even with these improvements, however, the method still 
requires cver 500 microseconds per variate in an assembler 
version. Thus, allowing for GAMA's substantial set-up time, 
it is faster to use the Johnk method for less than 250 
variates (k=0.75) to 1,350 variates (k=0.1), depending on 
k, but the method is too slow for large-scale simulation 
(100,000 or more variates). The method remains very 
atractive fcr use when time is not so critical, however, 
since it is very easy to program and requires minimal 


Storage. 
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EVEENGOIPULESCSPMULATION EXPERIMENT 


== et nt a а ње eee те ee ee ce Se es s w = > em AX E m MN 


ШІ DESCRIPTION OF METHOD 


The basic purpose of the simulation experiment was to 
mivestigate the distribution of Y for small values of J. 
J 


Different simulations were run for shape parameters X of .1, 


Et —5- and .75 corresponding to coefficients of variation 
2 
ЖИЕ ОР (0, 4, 2 and 1.33. In addition, a simulation using 


an implementation of the Marsaglia exponential generator was 


Bun Lor the case k=1.0. 


22 Simulation 


r — e a ше журты A а кәм mn 


Ten thousand samples of J gamma variates each were 
generated for J = 10,30,50,100 and the Р statistic computed 
for every sarple. The first four sample moments were found 
and then the entire sample of 10,000 rS vas ordered in 


order to estimate the quantiles of the distribution (see 
ENbSectron  IV.B.2). This procedure produced 27 statistics: 
four sample moment estimates and 23 quantile estimates. The 
entire procedure was then repeated ten times and the mean 
and variance of each of the 27 estimates found over the ten 
sanples of each. 

Obviously a great many gamma variates were required 
pos carry out the simulation, nearly 20 million for each 
value o£ k. The size of the Simulation graphically 


demcnstrates the need for speed in the gamma generator. 
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Even with the relatively fast performance of the 
Assembly-coded generator, the total IBM 360/67 central 
processor time needeä to run the Simulation was about 60 


minutes for eacn k value. 
2. FORTRAN IV Program 


Шек тат тсе Was Wprtten in FORTRAN IV; a 
program listing is included on pages 96-98. The program and 
subroutines required nearly 89,000 bytes of core at run 
time, 40,000 bytes of which were needed For the 
10,000-element array used for quantile estimation. А sample 
of output frcm the program for J-30 is included on pages 
na- 76. 


PER ESTIMATION OF PARAMETERS 


>; Jio r bt ion Moments 


-——À A— r m x w = ~ то с < = w ee ш оь 


Ti EY represents the i-th observation of Y in 
J; i J 


• 
$ 


the simulation, then the r-th moment m of У can be 
T J 


estimated by 


E E 2 
S = M 


in the simulation, estimates of the first four moments 


(r=1,2,3,4) were obtained. In addition, the sample variance 


EM tie coefficients of skewness and Laie oO Sis were 
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calculated from the moment estimates. 


2 Quantiles 


— oo P erum се o =з» 


а. Quantile Estimation Considerations 


ciao distribution F(x) 15 defined 
a 


as the solution to the equation 
Pex} =a (Оса <>"). 
а 


Quantiles are unique for most continuous distributions but 


are not, in general, for discrete distributions. In the 
simulation, quantiles of Y were estimated for a = .001, 
J 


EINEN ООО 20 eer, 050, —:1(.1).9, .950, .975, 
me, .990, .995, .998 amd. .999., 

There are two principal methods for  quantiie 
estimation in simulations: the order statistic estimator 
and the Robbins-Nonro stochastic approximation techniques. 
(For a more complete discussion, see Goodman, Lewis and 
Robbins [11].) The order statistic method was chcsen for 
this simulation because of the ease of implementation. Let 


the N observations Y m T. Y DIY be  crdered 
des Jg J;N J 


БӘ that Y ү єз ү ; then Y HEIC 
(9;1) (J; 2) (J; N) (2:1) 
called the i-th order statistic of the sample and an 


estimate of the a-quantile y of the distribution of Y is 
a J 


<! 
Н 


y е 
а (J; taN- ) 


2 





b. Sorting Random Data 
Io toSe-timdbe the extreme quautiles of 


the distribution of Y it is necessary to deal with large 
Jj 


Samples of Y. Considerations of bias (see the next 
J 


Subsection) also call for large samples. AS mentioned 
previously, a sample size of 10,000 was chosen for quantile 
estimaticn. 

БӘРЕП ch а large samples» can introduce 
considerable overhead into the simulation, however, unless 
the Sorting method is carefully chosen. There is a 
bewildering variety of known sorting methods for various 
purposes (Knuth [15] devotes an entire text to the subject). 
The particular technigue selected for the simulation is due 
to Singleton [27]; a FORTRAN implementation of the method is 
capable of ordering a sample of 10,000 normalized type REAL 
numbers in less than four seconds. In the simulation, 
therefore, less than five per cent of the execution time was 


spent in sorting the samples. 


3. Errors and Bias 

The parameter estimates may not correspond to the 
Exc parameter values for three reasons: statistical 
uxraability, computation errors and estimator bias. 
Statistical variability is an inescapable reflection of the 
КООШО that a simulation is not a deterministic process, as is 
the summing cf an infinite series. Thus, the probability 
that the precise expected value of an estimator will be 
observed is very lou and may be zero for a continuous 
distribution. The law of large numbers, however, guarantees 
that the estimator uill converge closer to its theoretical 
value as more replications of the simulation are performed. 


ато ӘГГОг 16 a reflection of digita) 
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Rn uter —round-off errors, that is, for every computer 
system there is a floating-point constant e such that when 


Ebr [el 
1.0 * d = 1.0. 


ШЕП спе ІНМ “3606 сустеп, е - 9.5310 . Thus, if very many 


data elements are to be summed (as in determining the moment 
estimators) there may be some loss of significance as the 
accumulated sum increases. The problem will be even more 
acute if very skewed data are to be added, as in the 


determination of Y анау _. 
1 28 


Round-off errors in > and > ‚ however, tend to 
J J 


cancel out when Y 1s calculated. Оп the other hand, in the 
J 


simulation the observed moment estimators m меге 
Б 


consistently belou their theoretical values, thus reflecting 
considerable round-off error. In fact, m was less than m 
r B 
17 observations out of 20 for r=1 and in 14 of 20 for 
r=2. Computing the moment estimators using double precision 
* a | * -16 
arithmetic (thus decreasing the value of e to 2.22 s 10 
NM ucrefore greatly reducing round-off errors) resulted in 
values of m more in agreement with theory. 
r 
Estimator bias 15 а theoretical property of the 
specific estimators used in a simulation. In this case, the 


Bere unbiased estimators of the m but 
E r 


ЕГУ - y Bla), 
LY 1 ү + sn 
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where N is the sample size (Goodman, Lewis and Robbins 


[11]. Thus y is a biased estimator. No attempt was made 
a 


Ко evaluate the "extent of the bias of y , but it is clear 
a 


that the simulation results may not be accurate for more 


than two or three decimal places. 


ЕЕ) SLMULATICN RESULTS 


1. Agreement with Theory 
For each simulation run, a t-statistic was computed 
using the sample mean and variance. Since the sample mean 


is the sum of 100,000 identically distributed random 


variables (the observations of Y ), it is safe to assume 
J 


Bat xts distribution iS very nearly normal and that 


t= (m, - E[Y J) / (S / VR) 


where 


_2 e 2 

Eu n (m — m ) 

п-т 2 1 

has a t distribution with n degrees of freedom. When n 15 
large (100,000 in this case) then the t distribution can be 
approximated by the standard normal distribution. The t 
statistics for the various runs are summarized in the 


following table: 
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J 10 30 50 100 


k 

(72 10 ӘС талог –0.5641 -0.9348 
ШЕ ІЗ II 51 $—2.40309 -1.9453 
0250 коп ОСОО 7-225376 1.2218 
0775 E2950. —2.6556 -1.2765 


1.00 о 006 -152623 МА 


Since the .01 quantile of the standard normal distribution 
is -2.326, it can be readily seen that the above results do 
not agree with theory. 

As mentioned above, the reason for this discrepancy 
Шо in accumulating rownd-off errors in calculating the 
moment estimators, The simulations for J=30 were thus rerun 
With the moment estimators computed in double precision. No 
change in the quantile estimates was observed, but the new t 


PratiStics vere 


k Ue Vc а 0.50 Ces 
estatistic = 0 2006 USO = 10500562 153176 


Thus, it is concluded that the simulation results agree with 


theory as far as the observed mcments are concerned. 


~ 


2: ПЕШИ по the Distribution of Y 


-—Ó — 


d 


тео ео the "actual values for У , the 
J 


simulation program produced normalized values for the 


quantile estimates 
E о) 715. 
à a 


where 
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|+ 
la 


and 


N) 


O = Var[ Y ] 
J 


E 


=a) = J+] 
12 Кә+1 


Mw norrnalized Statistics are tabulated here in crder to 


Show the asymptotic convergence of the distribution of Y to 
J 
a normal distribution; quantiles of the standard normal 


un bution are also listed. 





Alpha J-10 J= 30 4-50 15100 Normal 
Quantile 
50104 EA = 57515 20515 а О 
2002 = 2.1020 - 2.6114 -2.7360 -2.7978 = 22875 
„005 -2.1473 - 2.4244 22.49240 PO -2.996 
.010 -2.0788 SD O xcu -2.3166 22.320 
2020 -1,9445 -2.0244 -2.0907 22:090 229051 
E925 -1.8738 et) 759055 - 1.9672 = 12960 
„050 =1.6546 -1. 6654 - 16695 -1.6547 SION 
. 100 - 1.3444 03722 = У 0 54122905 2150757 
.200 -079230 -0. 8844 2:0 09897 028503 205082 
2500 20.5885 -0.5511 -0.5441 -0. 535! 0524 
.400 20522793 2072/06 -0. 2608 2022000 2022938 
2500 ООО -070051 №0923 020029 0 
.600 0.2802 0. 2687 0. 26095 0.2604 032933 
2700 0.5923 0.5555 059535 0.5382 0.524 
2300 2757257) 0° 56783 ЕЕ 2 022566 0. 882 
. 900 3306 Taa 1-3042 1:5 2958 12282 
950 1.6541 1.205953 1.00 2\ 1-6577 1.645 
2375 122754 1. 3228 9580 1: 9627 12:900 
2980 1723437 2.0036 2.0467 2. 0525 250 9H 
890 2a 2. 2905 2:2 2/9099 2+ 30916 2« 326 
25:05 2. 1498 2.4068 2. 5013 2.5484 2996 
os 2.1901 2.: 6027 227520 28258 22878 
2999 2.2009 DOO 2.9476 3.0085 3.090 
u E но 2/565 50.5 
С 22.035710 4. 3277 5.8914 8.7034 
Bot – 0. 00688 - 242118 moo 20999515 
Table I. Simulation results for kz0.10.  Normalized 


quantiles of the distribution of Y under the null 


hypothesis b=0 (no trend) are tabulated. 
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Alpha J=10 2530 2550 12100 Normal 
p Quantile 
.001 = 06865 oe ian О 2390909 АКО 
2002 -2.5614 ЕО =2.0066 = 2.0 2 = 2208108 
2005 =2 0792 = 55 5 -2.5519 22:0] = 576 
.010 — oe oU = -2.3405 5222705 
. 020 =2 0089 =2. 0596 ӘЛ «ШЕТІ? 522 2077 202052 
2025 -1.9340 ZN. MP 5 - 12960 
.050 = ti. Orem = 156725 ~a 6039 = 2300885 = 605 
„100 ЕВ uu о = 1229602 a 
.200 =a Gat 020063 USO 209950505 -03 94 2 
.300 -0.5511 -0. 5414 ИЗИ 2095999 - 0.525 
„100 -0.2643 - 0.2608 =03 2632 2057505 =07 203 
.500 0. 0053 -0.0016 -0.0070 QNS 0. 
„600 052890673 0. 2595 0,208 0. 2508 0.298 
2700 0.5609 0.5365 075311 0.5351 0.528 
. 800 0.8913 029639 O. 5 9 08529 0.314. 
. 900 769 1. 3024 1223075 1. 29%6 2222 
E950 12.6906 1:6636 26593 1.6576 1.615 
E975 19329 1220692 1550692 19776 1. 260 
749 2051091 2009099 2 09453 2.0682 2.051 
. 990 2.2295 22991 2.399 2.3398 2.326 
2552 2 29901 245438 2.5685 2.58459 2.576 
2:98 2.5484 2.196; 2 ОДЕЛО 2. 8870 Pas 2 
1799 2.6496 29108 3.0809 RONS 3.090 
u ЕЕ 15.5 2555 S ME 
O I5 153 2.906098 329276 СО 
But 120851 - 2758 = 220 zb. UNES 
Table II. Simulation results for k=0.25. Normalized 


quantiles of the distribution of Y under the null 


hypothesis b=0 (no trend) are tabulated. 
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Alpha О Jz30 2550 15100 Nornal 
B Quantile 
2:001 А 230597 ие =3. 0890 = 3. 090 
EOS 2708 25700 -2.9104 -2.9028 -2.878 
210805 e UIS -2. 5669 -2.6023 АО АЕС 
„010 -2.2861 — 72410 223116 59 о 
2020 =° ¿0 416 — 50051 220023 20925 = 25054 
2025 ee -1.9918 -1.9941 SIN = leo 
2050 - 199617 ЕОС ео - 1.0005 > 
2100 21.3020 = 13 100 Се ИР = 122986 OZ 
E200 ESO o> Ons ооо ООВ БОЕ 
2300 D -0.5460 722515 о! =01 528 
.400 EUM o0 ОС ОСО DEDIT IET Hon 2022953 
5200 OUS - 020051 00022 0: 0092 Oz 
2600 002665 02532 072596 05 2708 Ome oS 
. 100 0583 025375 0259,80 0. 5480 0.524 
=800 0. 8727 0. 8610 0. 8622 078679 Qu 2 
. 900 1. 3098 1.308 / 1.302 EZ 1,7202 
5550 16615 1.6212 1.6747 1.6861 1.645 
E975 1229576 129047 1.9851 220026 13260 
„980 2999 200855 3 220902 28090192 2.054 
.990 2.2750 288 2.9106 2.3790 22326 
83/95 2.23 223856 285955 2.6282 2576 
E98 2616926 238316 2.8901 2.9099 2.878 
.999 2. 90 98 3530650 3.0638 3. 1490 39090 

u eS те 2 50.5 

O 1. 1726 221639 2S 4.0421 

= st 0. 1989 = 326603 =2 85876 14 2289 
Table III. Simulation results for k=0.50. Normalized 


guantiles of the distribution of Y under the null 


hypothesis b=0 (no trend) are tabulated. 
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Alpha Jz10 3-30 2550 J=100 Normal 
EN Quantile 
.001 -2.9330 Ее 2o 1529 alt O 
(002 -2. 7684 =) =D 5 -2.9034 = 2.896 
2010 5 =2.5 166 270022 12954059 о д 
2010 =2 22957 =2. 3062 2022340 11223020 -2 9E 
2920 -2.0410 = 2.030 - 220684 2220875 - 2.051 
2025 — | 9544 1 | ее о 219907 = 1.200 
05:0 - 1.6617 -1.6686 oU o 21056678 bb 
.100 12099 И 24123063 IgE ID 
.200 SOTO =0.8556 ғ 25027 08497 -0. 842 
E00 =(). 5429 2085275 055391 Ole or io = 05521 
„400 -0.2628 =0-2%96 22000258 = 0 2475 о 3 
„500 0.0024 0. 0079 0.0003 0.0064 0 
- 600 0. 206898 0. 2658 02005 0.2631 07253 
„700 0.5410 0. 5380 0:55 355 022290 0.524 
.800 0.8707 0.8640 0.8574 0.8584 0.842 
.900 1203,53 29052 123030 1.3000 1299862 
29/550 120097 1. 687 1 1.6803 1.6708 AGUS 
"075 172598 1.9889 кооой к ОО 12960 
2980 20509 2.0905 2.0849 2.0894 2.051 
E990 22206 223905 2. 3639 2. 3518 2.320 
2995 25.02 246429 2.0022 225938 2 575 
998 2919333 2. 8948 229219 27 8904 2.575 
19 99 222501 ЗОЯ Э 3. 0848 >. 1195 3.090 
u 2472 13.5 2560 50,3 
c 0.9852 1.7855 NOONE | 83112 
EU -0.1440 —1. 8950 -2.6556 215295605 
Table IV. Simulation results for k=0.75. Normalized 


auanei l es or the distribution of Y under the mull 


hypothesis b=0 (no trend) are tabulated. 
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Alpha 4:10 1-30 4550 42100 Normal 
Quantile 

5001 -22 0558 LU IH ОО NA = j. 090 
1002 = m2 590 о 2 28018 
=005 -2. 5016 -2.5750 о =2, 596 
. 010 -2.2994 = 2.3304 E WS о 
020 - 2. 0442 SOT UO 2:20:35 xm TE 
8025 = "(129566 = 10655 ЕО 
„050 = .65 19 -1. 6464 246/499 - 1.645 
2100 zu. 2007 - 1.2846 = 2775 = 12872 
5200 -0. 8583 -0. 8440 -0.8348 250.902 
. 300 2027529 =). 5244 = 5172 20.524 
. 400 =0. 25%0 2:92 533 = 02409 2302 95 
500 = UNS -0.0004 0.0075 0 
.600 02562 052566 0.2607 0255 
“700 05286 0: 5298 05853 17 0.524 
. 800 0.8599 08422 0.8443 0.842 
9100 1: 2905 122875 1228943 1282 
#9 50 120535 TG 17 1. 6908 1.645 
B9 75 129521 Це 055 1736930 12960 
.980 22017 2.0588 2.0541 2. 054 
8990 291069 23152 231608 | 2:320 
E95 2:4933 2. 5500 DIOS 2509 16 
298 2007571 2998427 200796 22378 
2209 D 209 250438 SOUS 3.090 

u 9.5 15505 ZO 

O 0.8660 1 5906 220207 

К est 210525925 5C 2120523 
Table V. Simulation results for k-1.00.  Normalized 


quantiles of the distribution of Y under the null 
J 


hypothesis b-0 (no trend) 
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DEL ARGE-SCALE SIMULATION TOOLS 


In programming the simulation it became clear that a set 
of standard data-handling subroutines would have been most 
helpful for parameter estimation. The facilities of 
standard simulation languages (e.g., GPSS or CSMP) and data 
reduction programs (SNAP/IEDA or SPSS) are not suited to the 
amount or type of information generated öin large-scale 
Statistical simulations such as this one. Extensive 
research in the literature was thus required to find 
high-speed techniques for sorting, quantile estimaticn and 
Statistical testing and considerable effort was expended in 
programming the techniques. 

The sophistication of state-of-the-art methods in random 
number generation, for example, can result in substantial 
savings of computer time and core requirements but only at 
the expense of program conmplexity. Since the specific 
details of the method involved are usually of peripheral 
inportance to the simulation itself, removal of the 
programming burden will result in more resources becoming 
available for the main purpose of the simulation. 

following Lewis* discussion [24] of a large-scale 


Ewputer-aided statistical package, the following routines 
would have been of value in the Y simulation: 
J 
(1) А fully-documented quantile estimation 


subroutine using the maximum transformaticn Robbins-Monro 


Stochastic approximation technique. 


(2 One or more efficient  goodness-of-fit test 
subroutines: Ее ch) square test, Anderson-Darling 


Fest, Kolmogerov-Smirnov test. 


(3) A sorting routine equal in speed to ACH 
Шсегліпп 307 [27]. 
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TON 


(4) An efficient subroutine for accumulating moment 


estimators with minimal round-off error, 
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INNER LT ICABILITY OF RESULTS 


1. Gamma Generator 


In view of the wide variety of physical events 
which can ke successfully modeled empirically as gamma 
renewal processes, the availability of an exact, high-speed 
даппа variate generator should enable more extensive 
simulations to be conducted in many applications areas. ТЕ 
the shape parameter remains constant throughout a given 
 Ріісатіоп, pre-calculation of GAMA's tables and constants 
will save set-up time in succeeding runs; if not, then the 
trade-off between GAMA and the Johnk technique  (Subsection 
ив. 3) should be considered in selecting a generator. 
Nevertheless, it is clear that whenever a large number of 
gamma variates with shape parameter less than one is 
required, GAMA is superior to any other known method for 


generating them. 


2.  Iests for Trend 

Вере ОЕ О = section II.D indicate that the 
variance of the standard B test for the Poisson process 
needs to be "inflated" by the coefficient of variation 
Squared for use in a gamma renewal process. nnus ІК шау 
now be possible to accept the null hypothesis that a given 
process is without trend when prior to this result it was 
not. For example, Lewis and Shedler [20] in analyzing 


computer page exception data obtained B values which were 
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far too high to accept the null hypothesis (normalized 
eo ОГ 2499, 82607 aud —18.11 were observed; the .01 
ашалела 15 -2.53). Нонеуег, allowing for the extremely 
hdgh coefficients of variation involved (3.34, 3.27 and 
3.70), the null hypothesis can be accepted in two of the 
three cases. 

Since in practice the actual value of k is seldom 
known, the results obtained in Chapter II suggest that the 
estimated sample coefficient of variation could be used as a 
multiplier of the Poisson variance in applying the B test. 
The simulation results indicate that this procedure may lead 
to slightly higher acceptance levels in the gamma case when 
k is small, but the effect diminishes with increasing J. 

Without investigation of the relative power of the 
Ж test against other standard trend tests it may not be 
wise to apply it to an arbitrary process. It would appear 
most useful, in any case, when small values of k are 
involved, for this is precisely the situation which tie Б 
test does not handle well. 

As for the tabular results, Tables I - V are 
sufficiently accurate for routine application of the test 
based on B, especially since the level of the test is 
mebacrary. Possible inaccuracies in the tabular values 


might be significant, though, if the values were used to 
investigate the MUST of the E test against the 
И тог b#0. This has not been done. 

Ihe distributions of i are in all cases observed 


to be approximately symmetric and  underdispersed with 
Шо рест То the normal distribution. Convergence to the to 
Me normal distribution is rapid; the rate being slowest, as 
expected, for k=0.1. Even there, the normal approximation 
is adequate for J=50 if a level of no higher than 0.990 апа 


no lower than 0.010 is required. The approximate rates of 
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mommemcgence Of Y are such that rt appears possible to use 
J 


the normal approximation whenever J is greater than the 


following values: 


k 0.10 2225 050 O 75 1.00 
J 100 50 30 30 10 


B. AREAS FOR FURTHER STUDY 


1. | Gamma Generator 


The problem of a general purpose gamma generator (k 
> 1 as well as k < 1) should be investigated in order to 
extend the present work to wider applications areas. А 
straightforward way to do this is to follow Berman's method 
[28] of adding tk+ exponential variates to a gamma variate 
with shape parameter k ~ ‘К. This approach has two 
distinct draubacks in the present case: 

(1) When k - tk“ is small, GAMA requires much more 
Шор and generation time than when 0.5 = К - “Кі < 1.0; 
when К - tk < 0.05, GAMA cannot be used at all. This 
difficulty could. be overcome by using the fact that the 
square of a standard normal random variable has the gamma 
EE LIIbution ath k=0.5, Х=0.5; thus, if Z is a standard 
normal variate then Y - 22^ isa unit gamma variate with 
Shape parameter 0.5. ТЕ a Marsaglia normal generator is 
available, Y can be readily generated and added to a series 
of exponential variates to allow САМА to be called with 
Shape parameter in the optimum range. À thorough 
investigation of the time requirements for the various 
options is needed to find the fastest method in any given 


Case. 
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(2) When k is large, generation of tk} exponential 
variates may require excessive time. In this case, taking 
the logarithm of the product of tk4 uniform variates may be 
a more efficient way to proceed. | 

It may be possible to overcone bo t h these 
difficulties by extending the decomposition technique to 
values of k greater than one. In this case, however, the 
density СЕ Во опсес monotone and automatic 
generation of the sub-distributions and their probabilities 


Will be considerably more difficult. 
2. fests for Trend 


As mentioned previously, a more thorough 
investigation of the theoretical properties of the У test 
J 


is also indicated. Its relative power compared to other 
tests against the trend model should be determined; possible 
examples of such tests include a normal theory regression 
test for trend after a logarithmic transformation or the 
various non-parametric tests for trend that are essentially 
equivalent to Kendall's rank correlation test. Results 


Should also be obtained for other types of renewal process. 
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